Tensor network study of the Shastry-Sutherland model in zero magnetic field 
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We simulate the Shastry-Sutherland model in two dimensions by means of infinite projected entangled-pair 
states (iPEPS) - a variational tensor network method where the accuracy can be systematically controlled by 
the so-called bond dimension. Besides the well established dimer and Neel phase iPEPS confirms the presence 
of an intermediate phase with plaquette long-range order, and we determine its phase boundaries with high 
accuracy. The first order phase transition for J = 0.675(2) between dimer and plaquette phase is compatible 
with previous series expansion results. iPEPS predicts a weak first-order phase transition between plaquette and 
Neel phase occurring for J = 0.765(15). We do not find a stable intermediate columnar-dimer phase, even 
when we bias the state towards this order. 

PACS numbers: 75.40.Mg, 75.10.Jm, 75.10.Kt, 02.70.-c 



I. INTRODUCTION 

The Shastry-Sutherland model 1 (SSM) has been subject of 
many theoretical studies over the past three decades. Particu- 
larly after the discovery of the two-dimensional spin-gap ma- 
terial SrCu2(B03)2, 2,3 which can effectively be described by 
the SSM, 4 theoretical efforts to determine its phase diagram 
have been intensified. Since the model is frustrated, accu- 
rate Quantum Monte Carlo simulations are lacking due to the 
negative sign problem, and various other methods have led to 
conflicting conclusions so far. 

The SSM is given by the following Hamiltonian 



H = J^2s i -S j + J' Si 
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where Si are spin- 1/2 operators. The first sum goes over 
nearest-neighbors on a square lattice, and the second sum over 
next-nearest neighbor sites on orthogonal dimers following 
the pattern shown Fig. 1. In this work we set J' = 1 and study 
the phase diagram as a function of J. For J = the model 
reduces to a Hamiltonian of decoupled dimers with a ground 
state energy given by a product of singlets with an energy of 
—3/4 per dimer. This state remains the exact ground state also 
for finite J up to a certain value J c i. 1 The other limit J — > oo 
(or J' = 0) corresponds to the Heisenberg model, where the 
ground state exhibits Neel order. 

One of the first studies based on Schwinger boson mean- 
field theory predicted an intermediate helical phase between 
the dimer and the Neel ordered phase. 5 Other early works sug- 
gested a direct transition between the two phases without any 
intermediate phase. 4 ' 6 ' 7 

A plaquette phase as an intermediate phase was first found 
in the series expansion study by Koga and Kawakami. 8 They 
predicted a first order transition between dimer and plaquette 
phase for J c \ = 0.677(2), and a second order phase transi- 
tion between plaquette and Neel phase for J c2 = 0.86(1). 
This phase has been confirmed in the series expansion study 
in Ref. 9, where it was shown that the plaquette phase is adi- 
abatically connected to the ground state of the 1/5 -depleted 
square lattice model. A plaquette phase adjacent to the Neel 
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FIG. 1 : (Color online) The phase diagram of the Shastry-Sutherland 
model as a function of nearest-neighbor coupling J {J' — 1), ob- 
tained with iPEPS. The width of a bond is proportional to the mag- 
nitude of the bond energy, where full (dashed) lines correspond to 
negative (positive) energies. The arrows in the right panel illustrate 
the Neel order. In between the well established Dimer and Neel phase 
we find a phase with plaquette long-range order. 



phase has also been found in the theoretical study in Ref. 10 
based on a 1/N expansion. Further support for a plaquette 
state was provided by exact diagonalization and a combina- 
tion of dimer and quadrumer-boson methods, 11 where charac- 
teristic features of the intermediate plaquette phase have been 
found for 0.678 < J < 0.702. 

On the other hand, the series expansion study in Ref. 12 
challenged the prediction of a plaquette state as intermedi- 
ate phase, and proposed a columnar-dimer ordered state as 
another possible candidate. However, since none of the pro- 
posed non-magnetic states has an energy that is clearly lower 
than the Neel energy (in the relevant coupling regime) they 
were not able to make a final conclusion regarding the nature 
of the intermediate phase. 

In this work we show that the intermediate phase has in- 
deed plaquette long-range order, and we determine the phase 
boundaries with a higher accuracy than in previous stud- 
ies. Our results, summarized in Fig. 1, have been ob- 
tained by means of infinite projected entangled-pair states 
(iPEPS), which is a variational ansatz where the wave func- 
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tion is represented by a tensor network. 13,14 This approach 
can be seen as a two-dimensional generalization of matrix- 
product states (MPS) - the underlying variational ansatz of 
the famous density-matrix renormalization group (DMRG) 
method. 15 The accuracy of the ansatz can be systematically 
controlled by the so-called bond dimension D. A similar ap- 
proach has already been successfully used for the study of var- 
ious bosonic, frustrated spin- and fermionic models (see e.g. 
Refs. 16-22). We note that a finite PEPS has been used in a 
previous study of the SSM, 23 however, only with a small bond 
dimension D = 2, and no intermediate phase has been found. 

Besides the phase diagram of the SSM, one of the goals of 
this paper is to provide further benchmark data which demon- 
strate the performance and usefulness of iPEPS. We present 
results obtained with different simulation setups which are ex- 
plicitly biased towards certain orders. However, despite the 
bias, all simulation setups lead to consistent results in the large 
D limit. For example, we put a bias towards columnar dimer 
order, in which case for small D the dimer order is clearly re- 
produced. However, it vanishes for large D which shows that 
the dimer order is not stable. Finally, we present and test a 
scheme to treat next-nearest neighbor interactions in a more 
accurate way than in a previous study. 24 

The paper is organized as follows: In Sec. II we provide a 
brief introduction to the iPEPS method and explain the differ- 
ent simulation setups used in this work. In Sec. Ill we present 
our simulation results, first for values of J deep in the indi- 
vidual phases, followed by a detailed study of the phase tran- 
sitions. Finally, in Sec. IV we summarize our findings. In 
appendix A the scheme to treat next-nearest neighbor interac- 
tions in iPEPS is explained. 

II. METHOD 
A. Infinite projected entangled-pair states 

In this section we provide a short overview of iPEPS. For 
a more detailed introduction to iPEPS and tensor networks in 
general we refer to Refs. 14,25-27. 

The main idea of a tensor network ansatz is to represent 
(approximate) the coefficients c^^^.^of a wave function, 
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by a trace over a product of tensors. Here each index runs 
over the d local basis states of a lattice site. The most famous 
example are matrix product states (MPS) which form the class 
of variational states underlying the density-matrix renormal- 
ization group (DMRG) method. 15 In an MPS the coefficients 
are given by a trace over product of 3 -index tensors T\ r (with 
2-index tensors at the boundaries), as for example for a 6- site 
system 
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Thus, each coefficient cnwsiAwe * s gi yen by a product of ma- 
trices (with vectors at the open boundaries), hence the name 



FIG. 2: (Color online) Graphical representation of an infinite pro- 
jected entangled-pair state (iPEPS) made of a 4 x 2 unit cell of tensors 
(surrounded by thick dashed lines) which is periodically repeated. 
Each sphere corresponds to a rank-5 tensor and the lines (legs) at- 
tached to the sphere represent the indices of the tensor, as shown on 
the right hand side. 



matrix product state. Tensor networks are most conveniently 
represented graphically, as shown in Fig. 2(a) for this particu- 
lar example. Each tensor is represented by a shape with lines 
(legs) attached to it, which correspond to the indices of the 
tensor. A connection between two tensors implies a sum over 
the corresponding index, and an open leg of a tensor corre- 
sponds to the physical index for the local Hilbert space of a 
site. Each auxiliary index jh runs over D elements, which is 
called the bond dimension. Thus, D controls the size of the 
tensors (or matrices), i.e. the number of variational parame- 
ters of the ansatz. 

A projected entangled-pair state (PEPS) 13 is a natural gen- 
eralization of a matrix product state to two dimensions. In- 
stead of a three-index tensor, a five-index tensor T\ dru is in- 
troduced for each lattice site on a two-dimensional (square) 
lattice, where each tensor is connected with its four neighbor- 
ing tensors via the auxiliary indices I, d, r, u, each having a 
bond dimension D. Thus, the number of variational parame- 
ters per tensor is dD A . An infinite PEPS (iPEPS) is an ansatz 
for a wave function in the thermodynamic limit. 14 It is made 
of a unit cell of tensors which is periodically repeated on the 
infinite lattice, as depicted in Fig. 2(b). If the wave function 
is translational invariant, the same tensor can be used on each 
lattice site. If the state breaks translational symmetry, a larger 
unit cell may be required. 17 In practice, different unit cell sizes 
are tested to check, which size leads to the state with lowest 
variational energy. 

An iPEPS with D = 1 is nothing but a site-vectorized 
wave function (a product state), parametrized by vectors Ti 
on each site. With increasing D the iPEPS can represent more 
and more entangled states, with a scaling of the entanglement 
with block size which obeys the area law of the entanglement 
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entropy. 28 Or in other words, with increasing D the iPEPS 
can take into account more of the quantum fluctuations of the 
true ground state. These quantum fluctuations may select, e.g. 
one of infinitely many degenerate states in the classical D = 1 
case. Thus, iPEPS provides a way to systematically study 
a state as a function of D, where D controls the amount of 
quantum fluctuations (or entanglement) in the system. 

In order to obtain an approximate representation of the 
ground state for a given Hamiltonian, the tensors need to 
be optimized, i.e. the best variational parameters have to be 
found. In this work we do this optimization by performing an 
imaginary time evolution of an initial (e.g. random) iPEPS. 
The evolution operator exp(— f3H) is split into a product of 
two-body operators via a (second order) Trotter- Suzuki de- 
composition (see e.g. Ref. 27). Application of such a two- 
body operator to a bond in the iPEPS increases the dimension 
of the bond. Thus, for an efficient evolution, the correspond- 
ing bond needs to be truncated back to the original bond di- 
mension D after each time step. There are different schemes 
to perform this truncation. In the present work we use the so- 
called full update, which is explained in detail in Ref. 27 for 
Hamiltonians with nearest-neighbor interactions. In this trun- 
cation scheme, the full wave function is taken into account 
to determine the relevant subspace, in contrast to the simple 
update, 21,29,30 which is computationally cheaper, but less ac- 
curate. In Ref. 24 a simple update scheme for Hamiltonians 
with next-nearest neighbor interactions was presented. In ap- 
pendix A we introduce a more accurate full update scheme to 
treat next-nearest neighbor interactions. 

Once we have obtained an approximate representation of 
the ground state in form of an iPEPS we can compute expec- 
tation values of observables O. To do this we have to contract 
(evaluate) the tensor network representing (\I/|0|\I/). In an 
MPS this can be done in an exact way by performing a se- 
quence of pairwise multiplications of tensors. However, in 
iPEPS (and PEPS) the contraction cannot be done exactly (in 
polynomial time), and there exist different schemes to per- 
form the contraction approximately. Here we use the corner- 
transfer matrix method 31 ' 32 generalized to arbitrary unit cell 
sizes. 17 The accuracy of the approximate contraction is con- 
trolled by the "boundary" dimension which we typically 
choose between 20 up to several hundreds, depending on D. 
For the SSM, this is sufficiently large so that the error due to 
X is small (compared to the effect of the finite D). 

B. Order parameters 

Since iPEPS describes a wave function in the thermody- 
namic limit symmetries such as SU(2) spin rotation or lat- 
tice symmetries may be spontaneously broken. To check for 
SU(2) breaking we compute the local magnetic moment on 
each site, 

m = ^VW 2 + K> 2 + K>2. (4) 

A finite m implies that the SU(2) symmetry is spontaneously 
broken. 
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FIG. 3: (Color online) The four different simulation setups used to 
simulate the Shastry-Sutherland model with a square lattice iPEPS. 
Dark shapes correspond to tensors and dark lines to auxiliary bonds 
between tensors (the open, physical index of each tensor is omitted). 
Filled circles correspond to physical lattice sites and interactions be- 
tween the physical sites are given by thick shaded lines. In the dif- 
ferent simulation setups different sites are blocked together. A tensor 
is used for each of these block of sites: (a) Setup O: one tensor per 
physical lattice site (no blocking), (b) setup P: one tensor per four 
sites on a plaquette, (c) setup D: one tensor per orthogonal dimer 
(two sites), (d) setup C: one tensor per dimer (two sites), arranged in 
columnar order. 



To check for lattice symmetry breaking we measure the lo- 
cal bond energies on each bond in the unit cell, and then 
compute the following order parameter: 



AE = max(^) - min(E b ). (5) 



If AE is finite, then the state breaks lattice symmetry. By 
plotting the distribution of the different E b in the unit cell, as 
e.g. in Fig. 1, we can distinguish between plaquette-, dimer- 
or other possible orders. 

We note that broken symmetries can also be an artifact of 
a too small D in iPEPS. It is therefore always important to 
check, how an order parameter behaves as a function of D. 
If it tends to a finite value in the D — >> oo limit, then there 
is a true, physical symmetry breaking (long-range order). If, 
however, an order parameter is strongly suppressed with in- 
creasing D and extrapolates to zero in the infinite D limit, 
then this indicates that the exact ground state does not break 
the symmetry. 
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C. Simulation setups 

An iPEPS is a very general ansatz in the sense that the 
Hamiltonian is essentially the only physical input of a simula- 
tion (besides the choice of the unit cell size). The method is to 
a large extent unbiased, except that simulations for low bond 
dimensions are biased towards low-entanglement (mean-field) 
solutions, but this bias disappears for sufficiently large D. 
Nevertheless, in some cases it can be useful to bias the state 
towards a certain order, particularly in order to rule out a par- 
ticular order. For example, if we bias the states towards a 
dimer state, we are likely to observe a dimer state for small 
D. But if this order vanishes with increasing D, then this in- 
dicates that the order is unstable. 

We can put a bias by choosing different bond dimensions 
on the individual bonds in the iPEPS. For example, taking 
a much larger bond dimension between sites on a plaquette 
leads to a bias towards a plaquette state. Instead of using dif- 
ferent bond dimensions we can also block a certain number 
of sites together, and use one tensor for each block of sites. 
This effectively corresponds to having an infinite D between 
the sites within a block, i.e. all correlations between the sites 
in a block are taken into account exactly. We note that, no 
matter how we block the sites, in the large D limit one should 
always recover the same physics. Thus, trying different setups 
provides a way to check the robustness of a result. The same 
idea was used in Ref. 33, where we explained in detail how to 
construct the Hamiltonian for the block sites. 

In Fig. 3 the four simulation setups used in this work are 
presented. Setup O is corresponds to the "original" lattice 
with one tensor per physical lattice site. The horizontal and 
vertical shaded lines correspond to nearest-neighbor interac- 
tions, and the diagonal shaded lines to next-nearest neighbor 
interactions. Since there is no blocking and all bond dimen- 
sions are the same, there is essentially no bias in this setup. 

Setup P is biased towards a plaquette order by blocking four 
sites on a plaquette together, leading to a local Hilbert space of 
2 4 = 16 for each block site. One can see that there are only in- 
teractions connecting nearest-neighbor tensors (dark rounded 
squares). 

In setup D we use a tensor for two sites on each of the or- 
thogonal dimers, which leads to a square lattice iPEPS tilted 
by 45 degrees, with a local dimension of 2 2 = 4 for each 
block site. This setup is biased towards the orthogonal dimer 
state. We note that also in this case only nearest-neighbor in- 
teractions between the block sites appear. 

Finally, in setup C we also group two sites on dimers to- 
gether, but in a columnar-dimer arrangement, which puts a 
bias towards columnar-dimer order. The diagonal interactions 
oriented from bottom left to top right appear between nearest- 
neighbor tensors, whereas the other diagonal interactions con- 
nect next-nearest neighbor tensors. 

In setups C and O we make use of the next-nearest neighbor 
scheme introduced in appendix A, whereas in the setups P and 
D only a nearest-neighbor scheme for the interactions between 
the block sites is needed. 



D. Technical remarks on the iPEPS simulations 

There are different possibilities to choose an initial iPEPS 
for the imaginary time evolution. In many cases, an iPEPS 
obtained with the simple update provides a good initial state 
for simulations with the full update (cf. Ref 27). Typically, 
for small bond dimensions we evolve several random states 
and check which one has the lowest variational energy. This 
state is then used as an initial state for simulations at larger D. 
Close to a phase transition it can be useful in some cases to ini- 
tialize the iPEPS with a state deep inside a phase (away from 
the phase transition). For example, for the plaquette phase 
with setup D, we obtained particularly good results from an 
initial iPEPS which was obtained from simulations of the 
SSM with a plaquette bias (stronger interactions on plaque- 
ttes). 

In the imaginary time evolution the energy decreases as a 
function of /?. In some cases we observed a slight increase 
of the energy for large /3, probably because the evolution is 
only performed in an approximate way. In these cases we 
take the state with lowest variational energy, for the fj where 
the energy is minimal. 

Our simulation results are all obtained with a 2 x 2 unit 
cell. We tested larger unit cell sizes and did not find signs of 
another low-energy state which requires a larger unit cell. 

To improve the efficiency of our simulations we used ten- 
sors with Ij q symmetry (a subgroup of SU(2)). The tensors 
then acquire a block structure, similarly to a block-diagonal 
matrix. Details on the implementation of global abelian sym- 
metries can be found in Refs. 34,35. 

For the corner-transfer matrix scheme we adopted the 
method explained in Refs. 17,32, but instead of computing 
isometries based on a singular value decomposition to ab- 
sorb a column (or a row) of tensors into the boundary (see 
Refs. 17,32 for details), we use the projector introduced in 
Refs. 20,36. 



III. SIMULATION RESULTS 

Our main results are summarized in the phase diagram 
in Fig. 1: iPEPS predicts an intermediate plaquette phase 
between the dimer phase and the Neel phase in the range 
0.675(2) < J < 0.765(15). In the following we first dis- 
cuss the properties of the individual phases, and then provide 
a detailed study of the phase transitions. 



A. Dimer phase 

The dimer phase consists of a product of exact singlets 
along the diagonal bonds in the lattice, i.e. for each value J in 
the dimer phase the ground state is the same, with an energy 
of —3/4 per dimer, or E s = —3/8 per lattice site. 

With setup D this state has a trivial representation with 
iPEPS, i.e. it can be represented with D = 1 because the state 
simply corresponds to a product state of the dimers. For the 



5 



-0.66 
-0.662 
-0.664 
-0.666 
-0.668 

-0.67 



(a) 




o setup O / 
-□- setup P / 
— > — setup D I o 

* - setup C / 

MC •/ 

^ ■ 






0.1 


0.2 0.3 0.4 










1/D 




10" 1 










(c) 






i 


10" 2 






x ' - ■■ 




10" 3 










10" 4 




a 






10" 5 

( 










) 


0.1 


0.2 0.3 0.4 

1/D 





0.42 
0.4 

0.38 
E 0.36 

0.34 

0.32 
0.3 




0.1 0.2 0.3 0.4 

1/D 




-0.572 
-0.573 
-0.574 
. -0.575 
-0.576 
-0.577 
-0.578 



(a) 



o setup O 
-Q- setup P 
setup D 
setup C 



0.05 0.1 0.15 0.2 
1/D 



0.1 0.2 0.3 0.4 

1/D 



0.1 
0.08 
0.06 
0.04 
0.02 




(c) 










B var-~~'~~ 




□V 











0.3 
0.25 

0.2 
0.15 

0.1 
0.05 




(b) 



0.1 0.15 

1/D 



0.05 0.1 0.15 0.2 
1/D 




FIG. 4: (Color online) Benchmark results for the Heisenberg model 
( J' — 0, J — 1) for the different simulation setups, compared with 
state-of-the-art QMC simulations, (a) Variational energy as a func- 
tion of the inverse bond dimension, (b) Local ordered moment m as 
a function of inverse bond dimension. The linear extrapolations are 
a guide to the eye. (c) Deviation of the iPEPS energy from the QMC 
result, (c) (d) The difference in bond energies AE vanishes in the 
large D limit for all setups, which shows that lattice symmetries are 
unbroken. 



other simulation setups D = 2 is required, because a singlet 
cannot be written as a product state of two sites. 



B. Neel phase 

The Neel phase in the Shastry- Sutherland model is adiabat- 
ically connected to the Neel phase of the Heisenberg model, 
which corresponds to the limit J — »> oo (or J' = 0). In this 
limit the model is no longer frustrated and can therefore be 
solved by Quantum Monte Carlo (QMC). In Fig. 4 we com- 
pare the iPEPS results obtained with the different simulation 
setups with state-of-the-art QMC data from Refs. 37,38. Fig- 
ures 4(a) and (c) show how the variational energy gets im- 
proved with increasing bond dimension D. Our best vari- 
ational energy E^ =9 = —0.66939, obtained with setup D, 
agrees up to four digits with the extrapolated value from 
QMC, E QMC = -0.669437(5). 

In all setups the local magnetic moment m decreases with 
increasing D, as shown in Fig. 4(b), and approaches the QMC 
value, raQMc = 0.30 743(1). 38 It is not known how m depends 
on the bond dimension D which makes accurate extrapola- 
tions to the infinite D limit difficult. Empirically, a linear ex- 
trapolation for the largest few values of D yields a reasonable 
estimate of the value in the infinite D limit, with a relative 
error of the order of a few percents. This is of course much 
less accurate than QMC, however, we can obtain results with a 
similar accuracy also for the cases with finite J' , where QMC 



FIG. 5: (Color online) Results for the Shastry-Sutherland model in 
the Neel phase, (a)-(c) Energy, local order moment m, and differ- 
ence in bond energies AE as a function of inverse bond dimension, 
obtained for J = 1 with the four simulation setups. The linear ex- 
trapolations are a guide to the eye. (d) Local ordered moment as a 
function of J in the Neel phase, obtained with setup D. 



suffers from the negative sign problem. What really matters in 
the present study is that we can clearly distinguish between a 
finite and a vanishing order parameter to identify the different 
phases, which is clearly feasible with the current accuracy. 

Setup P (setup C) is biased towards a plaquette (columnar 
dimer) state, and this is why for small D one finds a finite 
value of the plaquette (columnar dimer) order parameter as 
shown in Fig. 4(d). However, the plaquette (dimer) order is 
strongly suppressed with increasing D and vanishes for large 
D. Thus, even though individual setups exhibit a bias for 
small D, eventually they all become exact in the large D limit. 

For finite J' the system is frustrated because the spins on a 
diagonal bond would like to be antiparallel to each other (in- 
stead of parallel as in the Neel phase). This competition of 
interaction leads to a suppression of the local magnetic mo- 
ment m. For example, for J = 1 (Fig. 5) the local moment 
is reduced torn = 0.21(1), and it is further suppressed with 
decreasing J. The difference in bond energies AE also van- 
ishes in this case in the large D limit, indicating the absence 
of translational symmetry breaking. 



C. Plaquette phase 

In Fig. 6 we present simulation results for J = 0.7, where 
all setups consistently predict a plaquette phase for large 
bond dimension. The lowest variational energy obtained is 
£f =1 ° = -0.3881 with setup D. 

The difference in bond energies AE in Fig. 6(c) clearly 
remains finite in all simulation setups, and we find that the 
strong bonds form plaquettes as illustrated in Fig. 1 for large 
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FIG. 6: (Color online) Results for the Shastry- Sutherland model in 
the plaquette phase, (a)-(c) Energy, local order moment m, and dif- 
ference in bond energies AE as a function of inverse bond dimen- 
sion, obtained for J = 0.7 with the four simulation setups. The lin- 
ear extrapolations are a guide to the eye. (d) Plaquette order param- 
eter as a function of J in the plaquette phase obtained with setup D. 



D. The dependence of AE on D varies from one setup to an- 
other, nevertheless, a rather large value between 0.1 and 0.15 
seems to be compatible with all simulation setups. 

The local magnetic moment in Fig. 6(b) is finite for small 
D, which would suggest coexistence of plaquette and Neel or- 
der. However, m is clearly suppressed with increasing D and 
it is very likely to vanish in all simulation setups in the large 
D limit, which shows that the SU(2) spin rotation symmetry 
is not broken in the plaquette phase, as expected. 

Finally, in Fig. 6(d) the plaquette order parameter as a func- 
tion of J for different bond dimensions D is shown. Its mag- 
nitude tends to decrease with increasing J but remains finite 
in the whole plaquette phase. 



D. Absence of a columnar-dimer phase 

In Ref. 12 a columnar-dimer state was predicted as one of 
the possible candidates for an intermediate phase from series 
expansion calculations. With iPEPS we can also reproduce 
such a dimer state with setup C which is biased towards this 
type of order. For D = 1 a product of vertical dimers is 
obtained in this case, as shown in the right inset in Fig. 7. 
Upon increasing D, however, this dimer order is strongly sup- 
pressed with increasing D, and plaquette order is enhanced, 
i.e. two parallel horizontal bonds connecting two dimers be- 
come stronger, as shown in the left inset in Fig. 7. We define 
the following columnar-dimer order parameter 
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FIG. 7: (Color online) Columnar-dimer order parameter as a function 
of inverse bond dimension for J — 0.7. The insets show the bond 
energies for simulation setup C, for D = 1 (right side) and D = 5 
(left side), where the width of a bond is proportional to the magnitude 
of its energy with full (dashed) lines corresponding to negative (posi- 
tive) energies. A finite dimer order parameter can be found for small 
D. With increasing D, however, the dimer order becomes weaker, 
and vanishes in the large D limit. 



with Eb x and E^y the bond energies in horizontal and vertical 
direction, respectively. This order parameter can be finite for 
small D, but even for setup C, it is seen to vanish in the infinite 
D limit, as shown in Fig. 7. Thus, the dimer state can appear 
as a low-entanglement solution, but eventually quantum fluc- 
tuations destroy the dimer order and give rise to a plaquette 
order instead. 



E. First order phase transition between the dimer phase and 
the plaquette phase 

In this section we determine the transition point J c \ be- 
tween the dimer phase and the plaquette phase. The transi- 
tion is clearly of first order, as can be seen in the jump in the 
plaquette order parameter at J c \ from zero to a finite value. 

To find the transition value precisely we can make use of 
the hysteresis effect in the vicinity of a first order phase transi- 
tion: When we initialize a state in the dimer phase and tune J 
to a value slightly above J c \ the state will remain in the dimer 
phase (since the state is metastable). Similarly, a state initial- 
ized in the plaquette phase will remain in that phase when de- 
creasing J to a value slightly below J c \ . The phase transition 
occurs where the energies of the two states cross. Since the 
energy in the dimer phase is independent of J, E s = —0.375, 
the transition occurs when the energy of the plaquette state 
crosses E s = —0.375. 

In Fig. 8, the energy of the plaquette state as a function of 
J is shown, obtained with setup D which yields lowest vari- 
ational energies. Since with increasing D the energy of the 
plaquette state is decreasing, the crossing point J® obtained 
for a certain D is an upper bound for J cl . For example, for 
D = 10 we find J% =1 ° = 0.6768. If we extrapolate the en- 
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FIG. 8: (Color online) Energy of the dimer state (horizontal dashed 
line) and the extrapolated energy of the plaquette state (full line). 
The first order phase transition occurs at J — 0.675(2) where the 
two energies cross. 



ergy linearly in 1/D using the three largest values of D, we 
obtain J^^°° = 0.6730, which we take as an estimate for 
the lower bound, since the energy is seen to converge faster 
than linear in 1/D. Thus, we find for the transition value 
J c i = 0.675(2), which is in agreement with the series ex- 
pansion result 0.677(2) from Ref. 8, and the prediction from 
a higher-order coupled cluster method (0.677). 39 



F. First-order phase transition between plaquette phase and 
Neel phase 



Locating the phase transition between plaquette phase and 
the Neel phase is the most challenging part of this work. We 
find that the transition is of (weak) first order. As discussed in 
the previous section the transition occurs where the energies 
of the two states cross. However, unlike in the transition be- 
tween the dimer and plaquette phase, the energies cross at a 
very small angle, which makes it difficult to very accurately 
determine the crossing point. 

In Fig. 9 we plot the energies of the two states in the vicin- 
ity of the transition as a function of 1 /D obtained with setup 
O. For J = 0.75 the plaquette state has a lower energy than 
the Neel state for large bond dimensions. By increasing J the 
energy of the Neel state is lowered with respect to the plaque- 
tte state, so that for J = 0.76 they become essentially equal 
at large bond dimensions, and finally for J = 0.78 the Neel 
state has clearly a lower energy than the plaquette state. Since 
the two curves have similar slopes we do not expect another 
crossing of the energies at larger D. 

A similar result is obtained with setup D, shown in Fig. 10, 
where the two energies become similar for J = 0.77 for large 
D. Thus, these results suggest a first order phase transition 
occurring for J c2 = 0.765(15) which is much smaller than 
the predicted value 0.86(1) from series expansion. 8 
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FIG. 9: (Color online) Energies of the plaquette state and the Neel 
state in the vicinity of the first order phase transition between the 
two states, obtained with setup O. For J = 0.76 the two states have 
a similar energy for large D. For J = 0.75 (J = 0.77) the plaquette 
state has a lower (higher) energy than the Neel state. 
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FIG. 10: (Color online) Energies of the plaquette state and the Neel 
state in the vicinity of the first order phase transition between the 
two states, obtained with setup D. For J = 0.77 the two states have 
a similar energy for large D. For J = 0.76 (J = 0.78) the plaquette 
state has a lower (higher) energy than the Neel state. 



IV. CONCLUSION 

In this work we showed that the intermediate phase in the 
Shastry- Sutherland model is a plaquette phase, in agreement 
with several previous studies based on other methods. 8 ' 9 ' 11 Us- 
ing state-of-the-art iPEPS simulations we have accurately de- 
termined the phase boundaries of the plaquette phase. By us- 
ing different simulation setups we biased the solution towards 
different ground states. But despite the bias, all setups consis- 
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tently predict an intermediate plaquette phase, and we found 
that the previously suggested columnar-dimer phase is unsta- 
ble. 

This work provides a good basis for future studies of the 
SSM in a finite magnetic field, a subject which is still attract- 
ing a lot of attention. Indeed, despite a considerable theoret- 
ical effort, 4 ' 40-43 the actual sequence of plateaus of the plain 
SSM has not been definitely established yet, at least close to 
the transition to the plaquette phase, and to get reliable in- 
formation on the magnetization curve in this parameter range 
is an important step towards the interpretation of the numer- 
ous and partially conflicting experimental results reported for 
SrCu 2 (B0 3 ) 2 . 2 ' 44 - 47 

Our study further demonstrates the usefulness of iPEPS as a 
tool to study strongly correlated systems in 2D. By simulating 
a model directly in the thermodynamic limit we can minimize 
possible boundary effects, which is an advantage over DMRG, 
where typically long cylinders with a certain width W are sim- 
ulated. For small widths DMRG is clearly more accurate than 
PEPS, however, since the number of relevant states to be kept 
in DMRG scales exponentially with W, DMRG becomes less 
accurate than PEPS for systems exceeding a certain width. 
For the Heisenberg model this seems to be the case around 
W ~ 10, 48 where the accuracy in the energy for W = 10 and 
m = 3000 states is comparable to the accuracy of an iPEPS 
with D = 9 for the system in the thermodynamic limit (for 
a PEPS with a finite width W = 10 and D = 9 the accu- 
racy would presumably be better). Thus, the two approaches 
can be seen as complementary: DMRG provides very accu- 
rate results up to a certain system size, and iPEPS approaches 
the problem from a different limit: Instead of varying the sys- 
tem size the only essential parameter is D, which controls the 
amount of entanglement in the system. While in many cases it 
is understood how to perform finite size extrapolations, little 
is known on how to extrapolate quantities in D. Obtaining a 
better understanding on how quantities depend on D will be 
important to determine order parameters more accurately in 
future studies. 

Finally, it is remarkable to note that in the above mentioned 
example for the Heisenberg model the total number of varia- 
tional parameters is roughly three orders of magnitude larger 
in DMRG (MPS) than in iPEPS, which shows that a PEPS 
offers a much more compact description of a 2D wave func- 
tion than an MPS. The challenge, however, remains to further 
improve the methods to optimize and contract an (i)PEPS, so 
that even larger bond dimensions and larger accuracies can 
be achieved. A lot of progress has been made in the last few 
years, however, we believe that a method which fully exploits 
the potential of the ansatz is yet to come. 

Note added. — After completion of this work, we became 
aware of a new study of the SSM based on the multi- scale 
entanglement renormalization (MERA) ansatz in Ref. 49, in 
which conclusions similar to ours have been reached for zero 
magnetic field. In particular, both methods locate the tran- 
sition between the plaquette and the Neel phase significantly 
below the estimate of Ref. 8. The critical value for the phase 
transition between the orthogonal dimer phase and the pla- 
quette phase reported in that paper, J = 0.687(3), is also 



not far from our result, but it lies between our D — 3 and 
D = 4 estimates, hence definitely above the upper bound 
J° =1 ° = 0.6768 provided by iPEPS. This indicates that the 
MERA result is closer to the mean-field solution than the 
iPEPS result. 
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Appendix A: Nearest-neighbor full update 

In this section we describe the full update used to per- 
form an imaginary time evolution with a Hamiltonian includ- 
ing next-nearest neighbor interactions. The description is in- 
tended for readers who are familiar with the basic notions of 
tensor network algorithms. For a general introduction we re- 
fer to Refs. 25,26. The idea of the following scheme is es- 
sentially the same as for the nearest-neighbor full update dis- 
cussed in Ref. 27: We apply a local time evolution operator 
U = exp(-rHi) to the iPEPS that represents a wave func- 
tion to obtain an evolved wave function (in imaginary 
time), 

|*')=#zl*> ( A1 ) 

This increases the bond dimension between the sites, on which 
Hi is acting on, from D to some D' = kD. The crucial step 
is now to truncate the enlarged bond indices from D' back to 
D, which yields an approximate wave function of 1^'). 
To do this we minimize the norm of the distance between the 
two wave functions, 

lll*'}-|*}l| 2 = + ( A2 ) 

with respect to the parameters of the new tensors in \Stf). 

As an example, lets consider the interactions on the trian- 
gle made of the sites A,B and C in Fig. 11(a). The SSM has 
three two-body terms on this triangle, which we denote by 
Habc- The corresponding imaginary time evolution operator 
Uabc — exp(— tHabc) can be exactly represented by a ma- 
trix product operator (MPO) with a certain bond dimension k, 
as shown in Fig. 11(b), with tensors X, Y, and Z acting on 
sites A, B, C, respectively. By multiplying these MPO ten- 
sors to the corresponding PEPS tensors, e.g. Y to B as in 
Fig. 11(c), we obtain a new PEPS tensor B' with increased 
bond dimension kD on the lower and the right leg. The aim is 
now to approximate tensor B' with another tensor B with all 
bond dimensions reduced to D. For this we make an ansatz for 
tensor B as shown in Fig. 1 1(d), made of the original tensor B 
and Y, and two tensors v and p of dimension Dx kxD. Sim- 
ilarly, we make an ansatz for A and C as shown in Figs. 1 1(e) 
and (f), respectively. 
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With these ingredients we can represent the wave functions 
|^), 1^'), as iPEPSs, where the new wave function de- 
pends on the parameters of the tensors u,v,p, and q. In order 
to obtain the optimal new iPEPS we need to minimize 

F(u,V,p,q) = (^(u,v,p, q )\^(u,v,p, q )}-^M(^(u,v,p, q )\^ f )) 

(A3) 

where we omitted the constant term \^') in (A2). The ten- 
sor network for the first term is shown in Fig. 1 1(g). As usual, 
the surrounding corner tensors C& and edge tensors Tj are 
obtained from the corner-transfer matrix method. 17 ' 31 ' 32 They 
take into account the infinite wave function (times its conju- 
gate) surrounding a 2 x 2 cell of tensors. The representation 
of the second term looks similar, except that the tensors u,v, 
p, q are replaced by identities (straight lines). 

As in the case of the nearest-neighbor update 50 we can min- 
imize this sum either iteratively, or e.g. with a conjugate gra- 



dient method. In this work we used an iterative scheme, where 
expression (A4) is first minimized with respect to u* with the 
other tensors being fixed, i.e. 

J-;F(u : v : p,q) = ^ [u*Mu - u*K] = (A4) 

where M is the tensor network in Fig. 11(g) with u and 
u* removed, and K is the tensor network representation of 
(^(u,v,p,q)\^ f ) with the tensor u* removed. Thus, the mini- 
mum can be found by solving the system of linear equations 

Mu = K, (A5) 

which yields a new tensor u. One then proceeds with tensor 
v, keeping tensors u, p, q fixed, and so on. This procedure is 
repeated until convergence is reached. 
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FIG. 11: (Color online) (a) We consider three sites A,B,C on a trian- 
gle in the Shastry-Sutherland model (cf. Fig. 3). (b) The imaginary 
time evolution operator on the triangle represented as a matrix prod- 
uct operator (MPO). (c) MPO tensor Y multiplied to the PEPS tensor 
B. This product can be represented by another PEPS tensor B' with 
enlarged bond dimensions kD in the lower and right leg. (d) Ansatz 
for the new tensor B to approximately represent B' . (e-f) Similar 
ansatz as in (d) for A and C. (g) Tensor network representation of 
(^|^) (see text). 



